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ABSTRACT 

In protoplanetary disks, the inner radial boundary between the MRI turbulent 
('active') and MRI quiescent ('dead') zones plays an important role in models of the 
disk evolution and in some planet formation scenarios. In reality, this boundary is not 
well-defined: thermal heating from the star in a passive disk yields a transition radius 
close to the star (< 0.1 au), whereas if the disk is already MRI active, it can self- 
consistently maintain the requisite temperatures out to a transition radius of roughly 
1 au. Moreover, the interface may not be static; it may be highly fluctuating or else 
unstable. In this paper, we study a reduced model of the dynamics of the active/dead 
zone interface that mimics several important aspects of a real disk system. We find 
that MRI-transition fronts propagate inward (a 'dead front' suppressing the MRI) 
if they are initially at the larger transition radius, or propagate outward (an 'active 
front' igniting the MRI) if starting from the smaller transition radius. In both cases, 
the front stalls at a well-defined intermediate radius, where it remains in a quasi-static 
equilibrium. We propose that it is this new, intermediate stalling radius that functions 
as the true boundary between the active and dead zones in protoplanetary disks. These 
dynamics are likely implicated in observations of variable accretion, such as FU Ori 
outbursts, as well as in those planet formation theories that require the accumulation 
of solid material at the dead/active interface. 

Key words: accretion, accretion discs — instabilities — MHD — turbulence - 
protoplanetary disks 



1 INTRODUCTION disk is fully turbulent, and beyond which the disk exhibits 

a characteristic layered structure for some range of radii. 
Protostellar disks, unlike the classic archetypes associated 

with dwarf novae and X-ray binaries, are enormous, cold, In general, it is unlikely that such a configuration sup- 

and poorly ionised. In fact, significant regions within these Ports steady accretion onto the protostar (Gammie 1996, 

disks are so weakly coupled to the magnetic field that MRI but see also Terquem 2008). In fact, observations of YSOs 

turbulence fails to develop (Balbus and Hawley 1991, Blaes show a rich assortment of time-variable accretion, the most 

and Balbus 1994, Gammie 1996, Igea and Glassgold 1999, striking of which are the outbursts associated with FU Ori 

Sano et al. 2000). Current models of protoplanetary disk and EX Lupi disks (Hartmann and Kenyon 1996, Herbig 

structure consist of a turbulent envelope of plasma (the 'ac- 2008). In addition, measurements of disk outflows reveal 

tive zone') encasing an extensive body of quiescent gas (the quasi-periodic variability on the order of years to hundreds 

'dead zone') (Gammie 1996, Armitage 2011). In particular, of years that may arise from fluctuating accretion in the in- 

these models posit a critical inner radius within which the ner disk, precisely at those radii where we might expect the 

inner dead/active zone boundary (Lopez-Martin et al. 2003, 

Raga et al. 2011, Agra-Amboage et al. 2011). 

Theoretical models that describe FU Ori behaviour at- 
* Email: hl278@cam.ac.uk tribute an important role to the dynamics of this critical 

f Email: steven.balbus@lra.ens.fr inner interface. In particular, many models interpret out- 
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burst events as the abrupt ignition of the dead-zone in tur- 
bulence, and hence the temporary dissolution of the criti- 
cal boundary (Gammie 1996, Armitage et al. 2002, Zhu et 
al. 2009, 2010). On the other hand, local numerical simula- 
tions of the MRI near marginality — the relevant regime at 
the interface — also exhibit strong variability in the form 
of violent intermittency (channel flows) (Miller and Stone 
2000, Sano and Inutsuka 2001, Latter et al. 2009, 2010) and 
oscillatory behaviour (Simon et al. 2011). The simulations 
demonstrate that the MRI will not smoothly 'die away' as 
one crosses into the dead zone through the transition layer. 
Despite these considerations, the interface is assumed to be 
effectively static in a number of planet formation scenarios 
(Varniere and Tagger 2006, Kretke et al. 2009, Dzyurkevich 
et al. 2010). The fate of these theories is unclear if indeed 
the interface is the site of the regular and violent dynamics 
suggested by observations and simulations. 

Our paper is concerned with the location, radial mor- 
phology, and dynamics of this important disk feature. We 
concentrate on an under-appreciated ambiguity in the loca- 
tion of the interface, one that can potentially have conse- 
quences for the observations and processes discussed above. 
Gas at the midplane is thermally ionised by the central star 
to sufficient levels for the MRI to operate, but only out to 
a critical radius 7?JJ ad very close to the star (< 0.1 au). Be- 
yond this radius lies a broader region in which stellar ra- 
diation fails but ionisation by turbulent heating succeeds. 
However, this requires there to be enough ionisation in the 
first place to get the MRI turbulence initiated. If the gas in 
this region starts off cold and poorly ionised it will remain 
so; conversely, if it begins hot and turbulent it can sustain 
this turbulent state via its own waste heat. In other words, 
the gas is 'bistable': it can fall into one of two quasi-steady 
states, active or dead. Finally, beyond a second critical ra- 
dius 7?* urb (~ 1 au) neither stellar nor turbulent heating is 
sufficient and the disk requires non-thermal sources to ionise 
the midplane. Usually this second critical radius is taken to 
be the actual inner edge of the dead zone. This need not 
be the case. The interface could in principle fall within the 
bistable region at smaller radii, and it need not be well- 
defined nor time-steady. 

In order to study this physics we develop a reduced 
MHD model that describes the interpenetration of the tur- 
bulence and the thermodynamics. It generalises the evolu- 
tion equations of Lesaffre et al. (2009), which couple the 
thermal energy with the turbulent energy, and installs a 
spatial (radial) degree of freedom. The gross disk dynam- 
ics are hence represented by two coupled reaction-diffusion 
equations. The main prediction of this model is that the 
dead/active interface, though well-defined, is not always sta- 
tionary. It tends to travel radially inward or outward de- 
pending on conditions at its present radial location. Ulti- 
mately the interface comes to a halt at a special radius, 
some 65% of the second critical radius i?* urb . Moreover the 
interface is relatively broad, with the temperature transi- 
tion layer extending over roughly 1 au. Numerical simula- 
tions should be able to check these predictions, while also 
assessing the prevalence and importance of turbulent inter- 
mittency. 

The organisation of this paper as follows. In Section 2, 



we develop the concept that gas in the inner zone of proto- 
stellar disks is bistable. In Section 3, we present a reduced 
model that displays the fundamental behaviour of the dy- 
namical system. This behaviour is more fully elucidated in 
Sections 4 and 5, in which the model is explored and its 
predictions presented. Next, the theory is applied, in Sec- 
tion 6, to what we argue is a more realistic model for the 
inner region of protostellar disks. Finally, in Section 7 wc 
summarize our findings and draw our conclusions. 



2 THE INNER REGIONS OF A 

PROTOPLANETARY DISK: A BISTABLE 

SYSTEM 

This section develops the idea that gas in the inner disk 
can fall into one of two states: a hot turbulent state, and a 
cooler quiescent state. We obtain a criterion for the insti- 
gation of MRI turbulence and show at which radii in the 
disk this is satisfied if we consider (a) radiation from the 
central star, in the absence of turbulence, and (b) thermal 
ionisation from turbulent dissipation, in the absence of stel- 
lar radiation. The latter radius can be significantly greater 
than the former, opening up a liminal region whose proper- 
ties we discuss. For clarity, we focus on gas at the midplane 
only, but our calculations can be generalised to include the 
vertical structure of the disk. 



2.1 Criterion for the onset of MRI 

The MRI only functions when there is adequate coupling 
between the gas and the magnetic field. An estimate of the 
critical ionisation fraction necessary for this coupling can 
be obtained within the bounds of resistive MHD. Hall ef- 
fects are of comparable importance in certain regions of the 
disk (Wardle 1999, Balbus and Terquem 2001, Wardle and 
Salmeron 2012), but for the sake of clarity they will be ne- 
glected; we expect no qualitative change in our calculation. 
The shortest MRI-active mode works on the resistive 
scale: l v = t)/va, where r\ is the Ohmic resistivity and va is 
the Alfven speed. This lengthscale cannot be larger than the 
disk scale height H . So we have marginal MRI when l v — H 
which gives 



C = v A H/ V = 1, 



(1) 



where C is the Lundquist number. To turn this into a crit- 
ical ionisation fraction we take H ss 0.1R and the following 
expression for resistivity 



ri = 234x- 1 T 1/2 cm 2 s- 1 , 



(2) 



where x e is electron fraction and T is temperature (Blaes 
and Balbus 1994). Combining these equations gives the crit- 
ical ionisation level: 



(ze)crit = 1-89 x W~ 14 /3 1/2 Rlv, 



(3) 



where j3 is the plasma beta and -Rau is disk radius in units 
of au (see also Balbus 2011). At 1 au and /3 = 10 2 - 10 3 we 
have the fiducial limit of x e ~ 10 -13 . 

Ionisation in a protoplanetary disk can be caused by 



non-thermal sources, such as cosmic rays, stellar X-rays, ra- 
dionuclides, and energetic stellar protons (Stepinski 1992, 
Gammie 1996, Igea and Glassgold 1999, Turner and Drake 
2009), in addition to direct heating from either the central 
star or turbulent dissipation. Each source gives rise to a 
spatial ionisation structure in the disk. In the case of ther- 
mal ionisation, this structure follows directly from the disk's 
temperature profile, and hence we can associate a criti- 
cal temperature Tmri with the critical ionisation fraction 
(xe)crit. The critical temperature may be computed with the 
aid of an appropriate form of Saha's equation, noting that 
the ionisation is controlled by the low ionisation-potential 
alkali metals (Pneuman and Mitchell 1965, Umebayashi and 
Nakano 1988). A straightforward calculation (e.g. Stone et 
al. 2000 or Balbus 2011) shows that Eq. (3) translates into 



Tmri = 800 - 1000 K. 



(4) 



Throughout the paper we take Tmri = 900 K as our refer- 
ence. 



2.2 Ionisation from turbulent heating 

MRI turbulence, if instigated, can thermally ionise the gas 
with some of its waste heat. The MRI can thus positively 
reinforce the unstable state from which it springs. In the 
complete absence of radiation, the MRI produces enough 
heat to keep the gas sufficiently ionised out to the critical 
radius i?' urb (Kretke et al. 2009, Balbus 2011). Essentially, 
the turbulence thermalises the free energy locked up in the 
background orbital shear. However, the available free energy 
declines with radius, and beyond 7?* urb there is not enough 
to keep things sufficiently hot if thermalised (given realistic 
disk opacities). 

The critical radius i?* urb can be calculated via an alpha 
disk model, which supplies a monotonically decreasing ra- 
dial temperature structure for the disk. The radius at which 
point T — Tmri corresponds to i?): urb . Of course, beyond 
i?* urb the alpha model is no longer consistent, because accre- 
tion via the MRI no longer functions in this region at every 
vertical level, but as a first approximation the approach is 
adequate. Kretke et al. (2009) summarised their numerical 
results by the following analytic estimate: 



Rl urb ^0.52 M% 9 Ml /3 



10- 



-1/5 



cm 2 g~ 



1/4 



(5) 
where M_7 is the mass accretion rate scaled by 10 7 solar 

masses per year, the stellar mass M» is scaled by solar mass, 
and k is the Rosseland mean opacity. With the positive cor- 
relation between M and M* , the more massive the star the 
greater _R); urb (Kretke et al. 2009). In summary, the critical 
radius associated with MRI turbulence falls between about 
0.1 and 3 au. 



2.3 Ionisation from stellar radiation 

Consider now a passive disk, in which there is no turbu- 
lence or accretion. It is irradiated by the central star and, 
potentially, cosmic rays. Though non-thermal sources, such 



as stellar X-rays and cosmic rays, are dominant in the sur- 
face layers of such a disk at small radii (and throughout the 
disk at large radii) , midplane gas at R < 1 au is well shielded 
and ionisation rates are vanishingly small. As a consequence, 
the associated x e is tiny (given standard dust-grain chem- 
istry) and the MRI is unable to work (Igea and Glassgold 
1999, Sano et al. 2000, Ilgner and Nelson 2006, Wardle 2007, 
Turner and Sano 2008, Turner and Drake 2009). More im- 
portant to midplane gas at these radii is stellar thermal ra- 
diation. 

A number of sophisticated radiative transfer studies 
have calculated the temperature structure of a passive disk 
irradiated by a realistic stellar source (see Pinte et al. 2009 
for references). These calculations omit turbulent transport 
in the active surface layers, but in any case the simulations 
of Hirose and Turner (2011) show that vertical turbulent 
transport of heat is negligible. The disk is usually truncated 
at a defined inner radius and the midplane temperature cal- 
culated as a function of distance from this inner edge. As 
discussed by Dullemond and Monnier (2010), at first the 
temperature structure is controlled by the dust sublimation 
threshold. However, once the sublimation temperature is 
reached, and dust can survive, T drops rapidly with radius. 
The temperature profiles calculated by Pinte et al. (2009) 
suggest that the midplane gas falls below Tmri when it is 
more than 0.01 au from the inner edge. This result is remark- 
ably robust across a range of normal optical depths (10— 10 6 ) 
and different extant codes. This finding is reinforced by the 
more detailed calculation of Woitke et al. (2009) which in- 
cludes a self-consistent account of the disk vertical structure 
and a great deal more radiative physics. These calculations 
imply that R 1 ^ may plausibly lie between 0.01 and 0.1 au, if 
the disk indeed has an inner edge. But if the disk physically 
connects to the central star, then R^ may be less. Given 
these uncertainties, we set a (crude) upper bound on i?JJ ad 
to be 0.1 au. 



2.4 The bistable region 

The ordering R T c ad < i?* urb invites us to split the inner re- 
gion of a protostellar disk into three zones as sketched out 
in Fig. 1. At very small radii R < R^, the disk will always 
be sufficiently irradiated by the star to ionisation levels that 
sustain the MRI. Similarly, the surface layers will also be 
unstable. Heat dissipated from turbulence is not required 
in these regions. At larger radii R > -Rc Urb , however, the 
miplane gas will always be too cold and poorly ionised to 
sustain the MRI. On one hand, it is adequately shielded from 
the star's radiation field. On the other, even if we managed 
somehow to kickstart the MRI at these radii, it would even- 
tually die out because its waste heat is insufficient to sustain 
the necessary temperatures. 

This picture leaves an ambivalent region i?J! ad < R < 
-R' urb sandwiched between the other two. Here the gas can be 
either quiescent or turbulent. It is bistable. If the gas starts 
off cold and quiescent it will remain cold and quiescent: the 
star's radiation cannot instigate the MRI in this region on 
its own. But if the gas starts off hot and turbulent it will 
remain hot and turbulent: the MRI will dissipate heat and 
this heat will be sufficient to maintain temperatures above 



•0 




activG 



Figure 1. Cartoon of a three zone disk. The very inner radii 
R < -R£ ad Ri 0.1 au (and also surface layers) are sufficiently ionised 
by the star to always support MRI. The outer radii R > ij* urb w 1 
au are always too cold to support the MRI — the turbulence 
cannot heat the gas sufficiently in this region and it is well shielded 
from the star. The middle region, sandwiched between the other 
two could go either way — turbulent or quiescent. It is 'bistable'. 
Note that at sufficiently large radii R > 10 au (not included in 
the cartoon) the disk may again be fully MRI unstable. 



the critical limit. The gas can only transfer from one state 
to the other via a sufficiently large (nonlinear) perturbation. 
This is a non-trivial point with potentially important 
ramifications. It leaves ambiguous the actual location of the 
dead/active zone boundary, and moreover implies that there 
may not be a well-defined interface at all. For instance, the 
gas here may oscillate between the two limits, or exhibit 
a range of interesting dynamical phenomenon, by analogy 
with other bistable systems (Zeldovich and McNeill 1985, 
Murray 2002). More generally, nonsteady behaviour in this 
region could explain some of the low amplitude accretion 
variability observed, and will probably be implicated in the 
dynamics of outburst events, if they are indeed caused by 
the aperiodic dissolution of the dead- zone (Gammie 1996). 



3 A REDUCED MODEL FOR THE 

TURBULENT/THERMAL DYNAMICS 

Our goal is to describe the dynamics in the bistable region 
of the disk, in particular the interplay between the ther- 
mal/ionisation dynamics and the turbulent dynamics. To 
that end we construct a mean-field model that captures this 
qualitative behaviour correctly. Our reduced model is rather 
crude, especially in its treatment of the turbulence, yet it il- 
lustrates in a clean way some basic ideas that should govern 
more advanced models and simulations. 

We the define the turbulent amplitude of the disk gas 
to be 



A=I(pv 2 + b 2 /(8^)>, 



(6) 



where v and b are fluctuations in velocity and magnetic 
field, p is density, and the angle brackets denote an az- 
imuthal and vertical average over the disk, as well as an av- 
erage over short radial scales (see Appendix A). Both K and 
the gas temperature T are intertwined in the bistable zone: 
the magnitude of the turbulent amplitude will determine the 
temperature via turbulent dissipation, but the magnitude of 
the temperature will determine the ionisation fraction and 
consequently whether turbulence is present or not. These 



relationships make up the 'reaction terms' in our dynamical 
model, and have been investigated on their own previously 
by Balbus and Lesaffre (2008) and Lesaffre et al. (2009). In 
addition, we add the influence of turbulent diffusion, treating 
T and K as reactive scalars in a turbulent flow field which 
both transports and is reacted upon by T and K. 

In effect, we construct a turbulent closure model with a 
novel dynamical connection to the thermodynamics. Anal- 
ogous models have been used in various terrestrial contexts 
with some success, such as the k — e turbulence model (for 
example, Davidson 2000), the multi-scale fluid models (see 
Frisch 1995 for references), and in mean-field plasma turbu- 
lence (for example, Gruzinov and Diamond 1994). Of course, 
the closest cousin of the model we present is the alpha disk 
itself, as interpreted by Balbus and Papaloizou (1999), and 
its more complicated variants (Armitage et al. 2002, Wunsch 
et al. 2005, 2006, Zhu et al. 2009, 2010). The detailed deriva- 
tion of our governing equations we give in Appendix A; in 
the main body of the paper we motivate their main features 
heuristically. We begin their exposition with a summary of 
the closely related system in Lesaffre et al. (2009). 

3.1 Homogeneous dynamics 

Lesaffre et al. (2009) argue that the full equations of vis- 
cous and resistive MHD in a shearing box can be effectively 
modelled by the following simple dynamical system: 






A(T), 



(8) 



where K is the box-averaged turbulent kinetic and magnetic 
energy, and T is the box-averaged temperature. The evolu- 
tion of these quantities is controlled by three functions: Q, 
which summarises the growth and saturation of the MRI, F 
which represents turbulent heating, and A which represents 
radiative cooling. Though these functions can be written as 
complicated averages over the fluctuating quantities, they 
may be more simply modelled by physically motivated al- 
gebraic expressions. When quantitatively compared to full 
MHD simulations, this simple system does remarkably well 
in describing the principle dynamics (Lesaffre et al. 2009), 
which inspires confidence in its utility. 



3.1.1 Parametrising the MRI 

To describe the onset and saturation of the MRI via the 
term Q a simple Landau operator is invoked for a single 
mode amplitude, which yields 



QxsK-aK 1 



(9) 



where s is the MRI linear growth rate and a is a parameter 
associated with its saturation. The temperature influences 
the linear forcing term, the growth rate, via 



s = so[l-v(T)], 



(10) 



where so is the ideal MHD growth rate, and rj is the magnetic 
diffusivity scaled by a constant reference diffusivity. Lesaf- 
fre et al. use the Spitzer prescription for r\ but in weakly 



ionised protoplanetary disks one based on an appropriate 
Saha equation is required (outlined in Appendix A) . 



3.1.2 Turbulent heating 

The turbulence taps energy stored in the background differ- 
ential rotation. Turbulent stresses can transform this energy 
into fluctuations that tumble down a cascade (or pseudo- 
cascade) into the arms of the dissipation scales where it is 
degraded into heat. The rate of energy injection into the 
system is hence proportional to the quadratic correlation 
{vrV,/, — bjib < f,/(4:TTp)). Lesaffre et al. approximate this cor- 
relation by the turbulent magnitude K (itself a quadratic 
quantity in the fluctuations) and write 



equations (7) and (8). This augmented model would then 

resemble 



r = wK, 



(ii) 



where w measures the local shear rate. The greater the local 
shear the more energy the turbulence can extract. We hence 
expect w to be a decreasing function of radius, though in 
most applications we will leave it as a constant. 



3.1.3 Radiative cooling 

Finally, if we assume that the disk surfaces radiate like a 
blackbody, the cooling function A may take the following 
profile 



A = b(T 4 -T 4 q ), 



(12) 



where T cq is the temperature at radiative equilibrium and 
b is some free parameter associated with the opacity of the 
gas. It controls, in basic terms, the speed at which radiative 
equilibrium is enforced. The cooling time is hence ~ 1/(6T 3 ), 
which we expect to be generally much longer than an orbit 
in the bistable region. A small optical thickness corresponds 
to a large b and a large optical thickness corresponds to a 
small b. Consequently, b may be a function of R in general. 
Note that in Lesaffre et al. (2009) a simpler linear cooling 
was employed. 



3.2 Turbulent transport 

Equations (7)- (8) describe the average dynamics at a fixed 
radius in the disk, however our goal is to capture the disk 
evolution over a significant range of radii and when subject 
to turbulent transport. Certainly, MRI turbulence effectively 
moves angular momentum outward, but it will also transport 
heat, though perhaps less efficiently. The MRI will also tend 
to spread kinetic and magnetic energy: the more aggressive 
turbulent fluctuations in higher intensity regions will travel 
into less turbulently intense regions. The latter transport is 
analogous to the spreading of kinetic energy away from a re- 
gion of localised forcing and into unforced quiescent gas. The 
disordered motions decay as they spread, but these residual 
motions still have energy distributed over a range of scales 
which ultimately dissipates. 

The simplest way to model the transport of T and K 
is as a diffusive process. We may then permit T and K to 
vary with radius 7?, in addition to time, and subsequently 
erect Fickian diffusion (or 'eddy viscosity') terms in both 



dK _, 1 d ( n „ dK 



dT „ . 1 9 / nn 9T 

!n= r - A+ R8R [ RDt M 



(13) 
(14) 



in which we have introduced diffusion coefficients Dk and 
Dt- Note that the latter coefficient not only incorporates 
transport by turbulent eddies but also by radiation. Both 
Dk and Dt should be (increasing) functions of the turbu- 
lent intensity K, and Dt may depend on temperature as 
well. However, for simplicity, we will assume that they are 
constant. 



3.3 Dimensionless equations 

The number of free parameters in the system can be re- 
duced by choosing suitable dimensions. We scale time by 
1/so which is of order an orbit, and we scale turbulent am- 
plitude K and temperature T by reference values K eq and 
T cq . We denote by K eq the turbulent amplitude in the ab- 
sence of resistivity (recall that T oq is the temperature in the 
absence of turbulence). The quantities so, T eq , and to some 
extent K aq , may be associated with a prescribed disk radius 
lying within the bistable region, between 0.1 and 1 au. Next 
we set the diffusivitics 

Dt = L T so Dk = L K so 

where Lk and Lt are the associated diffusion lengths at the 
prescribed radius. We next scale the unit of length in our 
equations by Lt- Together these transformations yield the 
simpler set of reaction-diffusion equations: 



dK _„ r . 2 n 1 d f n dK\ 

1h= sK - k +Pr Rlm{ R lm) 



dT _ 



Here 



**-»^*l 



wK cq - _ bT cq 

SO T eq SO 



(15) 
(16) 



and we have defined the Prandtl type number P r — L K /L T . 
From Appendix A, we have a growth rate prescription ap- 
proximating the influence of the Saha ionisation law 



a = tanh [6(T - Tmri)] • 



(17) 



This means that near the temperature Tmri the ionisation 
fraction jumps rapidly to a value sufficient to instigate the 
MRI. If, in physical units, T oq ~ 250K (associated perhaps 
with R = 0.5 au), we have Tmri ~ 3.6. There are now 
three parameters w, b, and P r . For notational ease the over- 
lines over s, to, and b will be suppressed in the following. 
Finally, we adopt P r — 1, not only for simplicity but also 
because we have no numerical estimates for the relative ef- 
ficiency of turbulent K and T mixing. Thus Lt = Lk- The 
model offers no detailed quantitative constraints on the two 
remaining parameters w and b, but their meanings are rela- 
tively transparent, measuring shear rate and inverse optical 



thickness, respectively. Because the injection of energy via 
turbulence proceeds on a time scale longer than an orbit, 
we take w < 1. Similarly, the cooling time 1/(6T 3 ) must be 
longer than the orbital time. With T = Tmri, this yields 
b < 0.01. 



3.4 A further simplification: slavery 

It is possible to reduce the system (15)-(16) even further if 
we assume that (a) the time-scale of the MRI saturation is 
much faster than the thermal timescale, and (b) that the 
diffusion of K is slow compared to the diffusion of T, i.e. 
P r < 1. Neither assumption may be strictly true, but the 
simpler system is particularly convenient to analyse, allow- 
ing us an insight into how the more general system works. 

These assumptions mean that K is always in equilib- 
rium on the long thermal times that concern us. In other 
words the turbulent K dynamics are 'slaved' to the thermal 
T dynamics. For long times we have 



K 



if TsCTmri 
$ if T > Tmri ■ 



(18) 



Now we need only solve one evolution equation, which we 
write as 



dt + RdR\ or 



(19) 



where as before the cooling rate is A = b(T 4 — 1) and now 
the heating rate is the continuous piecewise function: 



r = 



!) 



if T ^ Tmri 



W s if T > Tmri 



(20) 



3.5 Caveats 

Before applying the model we should stress a few points re- 
garding some of its weaknesses, especially with respect to its 
simplistic treatment of the turbulence. As is made clear in 
Appendix A, the validity of the model rests on the soundness 
of averaging over the small-scale disordered motions of the 
MRI turbulence, in addition to the vertical thickness of the 
disk. In particular, it demands a clean separation between 
large scales of at least H, upon which global disk proper- 
ties manifest, and the short MRI turbulent scales. Stratified 
simulations of the MRI suggest this is only marginally true 
at best (Davis et al. 2010, for example), while near critical- 
ity the fastest growing mode varies on approximately H . We 
believe this problem will not derail the qualitative results. 

In addition, for the small scale averages to make sense, 
the fluctuations need to be well-behaved — they cannot be 
subject to intermittent and volatile outbursts or other crit- 
ical behaviour. Yet that is what we might expect from the 
marginal MRI! The best we can hope for is that this bad 
behaviour does not invalidate the predictions of the mean- 
field model — that, in fact, the model indeed does a decent 
job of describing the mean dynamics, even if the fluctuations 
around that mean are a little wild at times. 

One should also add that the eddy viscosity descrip- 
tion of turbulent flow has long been known to be at best a 



crude, and sometimes misleading, approximation. In reality 
the stress exerted by such flows are non-Newtonian, some- 
times exhibiting a nontrivial dependence on the shear rate, 
effects associated with their finite relaxation time, and even 
negative transport (see, for example, Frisch 1995, Davidson 
2000). That all said, we are interested only in the basic dy- 
namics, which we believe are independent of the particulars 
of the turbulence itself. 

A final issue is the vertical averaging: in doing this, we 
are discarding any information about how the disk thickness 
responds to local changes in either the turbulence or tem- 
perature. It is implicitly assumed that H varies passively 
with changes in T. For example, a transition between the 
hot MRI active region and a cold quiescent region will be 
characterised by a drop in H, meaning the inner turbulent 
radii of a protostellar disk may bulge. Any interesting sta- 
bility problems, or other dynamics, that might arise from 
this configuration cannot be captured by the model. 



4 ANALYSIS OF HOMOGENEOUS STATES 

4.1 Local states 

Our first task is to compute homogeneous steady states — 
i.e. equilibrium states with no radial structure, which de- 
scribe either dead or active zones. As in a shearing box, 
such equilibria we expect to hold locally in a protoplane- 
tary disk. Because the diffusion terms take no part in these 
calculations, we need only solve a coupled pair of ODEs. 
The steady-state equations to be solved are: 



sK-K- 1 



0, 



wK -b(T 4 - 1) = 0, 



(21) 



which arise directly from Eqs (15)-(16). The first equation 
matches the linear forcing by the MRI to its saturation 
through the quadratic term. The second equation matches 
the first turbulent heating term with the second radiative 
cooling term. 

This system admits the cold quiescent (dead) state 



K = 0, 



r = i. 



Turbulence is absent and the gas falls into radiative equilib- 
rium. The system also potentially admits turbulent states 
K > 0, with K = (b/w)(T 4 — 1) and with the temperature 
determined from the nonlinear algebraic equation 



ws = b(T 4 -1), 



(22) 



with s given by Eq. (17). Here the heating rate r = ws is 
on the left and the cooling rate A = b(T 4 — 1) on the right. 
Naturally, the linear growth rate s must be positive for this 
to work or else there is no turbulent heating at all. Note 
also that the equilibrium depends on a single parameter, 
the ratio b/w. 

Equation (22) must be solved numerically. But if we 
plot the heating and cooling as functions of T we can get a 
good idea of how the balances work. Examples are given in 
Fig. 2 for three representative values of b/w. The points at 
which the two curves cross correspond to turbulent /thermal 
equilibria. Note that, depending on the ratio b/w, we can 
achieve one solution (the cold quiescent state), two solu- 
tions, or three solutions. In the figure we plot these three 
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Figure 2. Heating and cooling as functions of temperature T 
for different ratios of b/w = 0.001, 0.0042, 0.01. Crossings cor- 
respond to thermal equilibria and are marked with black dots. 
These comprise the cold quiescent state T = 1, and up to two 
turbulent states: Tj (which we term the 'intermediate state') and 
Ta (which we call the 'active state'). The former is unstable and 
the latter is stable. 



different cases. The critical b/w value above which turbu- 
lence vanishes, and for which there is only the dead state, is 
approximately 0.0042. When b/w > 0.0042 the cooling rate 
is too efficient, and/or the amount of free energy insufficient, 
for the requisite MRI temperatures to be maintained. Natu- 
rally we associate the cold quiescent state with a dead zone, 
and the hotter of the two turbulent state with the active 
zone. The intermediate state we show later is unstable and 
will not be observed, though it is important in policing the 
boundary between the other equilibria. 

Another way to characterise the equilibria is to rewrite 
Eq. (22) as 

T=[l + (w/b)s] 1/4 . 

Because s is essentially a step-function in T, and the upper 
turbulent solution always occurs when resistivity has been 
virtually 'switched-off ', we can approximate the upper state 
T A by 



T A ^(l + w/b) 1/4 , 



K A 



(23) 



Therefore, the greater the parameter combination w/b the 
greater the active zone temperature Ta- This makes sense 
because w/b measures the amount of heat stored in the gas. 
Recall b represents the ability of the gas to hold heat (opti- 
cally thick gas takes smaller b) while w controls the amount 
of heat generated via turbulence (w is larger when more 
energy is available to be dissipated). Optically thicker gas 
at smaller radii is hotter, whereas optically thinner gas at 
larger radii is colder. 



4.2 Linear stability 

Next we briefly inspect the stability of the three possible 
steady states. We find that the quiescent state Kq — 0, 
Tq = 1 is always stable on account of its cooling law (heat- 
ing is absent around this state). The upper turbulent state 
Ka ~ 1, Ta ~ (1 + w/b) 1 ' 4 is also stable. This is easy to see 
from the heating and cooling functions in Fig. 2. If the tem- 
perature increases from the hotter active state the cooling 
function will be greater than the heating function: the ex- 
cess temperature will be removed and the system will return 
to equilibrium. On the other hand, if the temperature de- 
creases then heating will dominate cooling and the gas will 
also return to equilibrium. It is quite the opposite for the 
intermediate equilibrium: a small positive deviation in tem- 
perature means that heating will dominate cooling and so 
the gas will continue heating up. Similarly a small negative 
deviation means that cooling will dominate. 

This heuristic reasoning can be confirmed with a de- 
tailed linear stability analysis, which we present in Appendix 
B. The linear dispersion relation that ensues gives the fol- 
lowing stability criterion 



(24) 



dA dT 

dT > dT' 

which puts in mathematical terms the graphical argument 
presented earlier. Direct substitution of T = wK and A = 
b(T 4 — 1) shows that the quiescent state and the upper active 
state are always stable, while the intermediate state is not. 
It is, in fact, a saddle point with one stable and one unstable 
mode (manifold). 



4.3 Nonlinear homogeneous dynamics 

Lastly, we describe the nonlinear dynamics of the homoge- 
neous equations. This situation might mimic conditions in 
a localised patch, as approximated by a shearing box (cf. 
Lesaffre et al. 2009). The full dynamical system here is: 



dK 



= sK - K z 



dt 

§=wK-b(T 4 -l), 



(25) 
(26) 



two nonlinear ODEs in time, which we solve numerically 
using a Runge-Kutta method. In Fig. 3 a number of numer- 
ical solutions to this system are plotted alongside the vector 
'reaction' field R which controls the solution trajectories: 



R= [sK - K 2 , wK -b(T 4 - 1)]. 



(27) 



The three equilibrium points are represented by black spots. 

As is clear, trajectories eventually end up on one of the 
stable states. The speed at which the system moves to one 
or the other state is limited by the cooling rate ~ bT 3 . Most 
trajectories are first attracted to the unstable saddle point, 
guided along its stable manifold (the red dashed curve), but 
are finally repelled along the stable manifold (the dotted 
curve) . Which state they terminate upon depends on where 
in the phase space they begin and, more specifically, in which 
basin of attraction they initially fall. 

The phase space can be divided into two regions, the 
basin of attraction of the quiescent dead state, and that of 



also introduce an intermediate radial variable X so that 




Figure 3. Phase portrait of the homogeneous problem. Some so- 
lutions of Eqs (25) and (26) are plotted in blue along with the 
reaction vector field in black. The three equilibrium solutions are 
represented by black dots. In addition, the stable and unstable 
manifolds of the intermediate state are plotted as dashed and dot- 
ted lines respectively. Trajectories which begin above the dashed 
line (the stable manifold) end up at the turbulent state, and tra- 
jectories which begin below the dashed line end up in the dead 
state. The average time it takes to reach either state is limited by 
the cooling time. A trajectory settles on the hot turbulent state 
in some 10 orbits, whereas one approaches the quiescent state in 
over 100 orbits. Parameters are w = 1, b = 0.001. 



the active turbulent state. The boundary between the two is 
given by the stable manifold of the intermediate state (the 
dashed curve). As the parameters change, this boundary also 
changes, with the basins of attraction growing or shrinking 
relatively. In particular, as b/w decreases the intermediate 
state approaches the low state, thus diminishing the basin 
of attraction of the low state. 



5 MRI-FRONTS 

We have dealt with the homogeneous dynamics of our re- 
duced model, first detailing the equilibrium states the sys- 
tem admits, their linear stability, and the ensuing nonlinear 
dynamics. The analysis has been limited to gas at a par- 
ticular radius, allowing the gas to settle on either an active 
or dead state. But, as we vary the radius in the disk, the 
gas presumably transitions from one state to the other at 
some sort of interface. In this section we attempt to compute 
'fronts' connecting these two stable homogeneous equilibria. 
We hence must reinstate the spatial element of the nonlinear 
problem. 



X = R-R . 



The equations are now 
dK 



= sK-K z + 
at 



d 2 K 
dX 2 ' 



%- =wK -b(T 4 -I) 



8 2 T 
dX 2 ' 



(28) 
(29) 



Next we make the following ansatz: there exists a front with 
a radial profile that is steady and which moves radially at 
a constant speed c. This front connects the stable quiescent 
equilibrium state {Kq, Tq) and the stable turbulent state 
{Ka, Ta)- Because of radial disk structure, the assumptions 
of constant c and steady profiles only hold approximately 
in our small region near _Ro- However, they permit us to 
directly compute the fronts and to clearly examine their 
properties. We relax these assumptions later in numerical 
simulations. 

To make progress we move into a frame comoving with 
the front, and introduce the comoving spatial variable 



£ = X -ct. 



(30) 



In this frame the front profile is stationary and thus K — 
K{() and T = T(£). The partial differential equations (28)- 
(29) become the coupled ordinary differential equations 



d 2 K dK 



sK-K 2 



0. 



d 2 T dT „ w^4 

ie +c -i +wK - b{T 



1) = 0. 



(31) 



(32) 



The boundary conditions for these equations come from the 
requirement that the front must connect {Kq, Tq) with 
{Ka, Ta)- Formally at least, this may be written as 



K 
K 



Kq, 
K a , 



T 
T - 



Tq 
T a 



£ ->■ oo, 

£ — > —00. 



(33) 
(34) 



A more compact vector form for the equations (31) and (32) 
is 



R = 



(35) 



d z dz 

d^ + c di 

where we haved introduced the solution vector z = [K, T]. 
Recall that R is the reaction field introduced in (27). It 
corresponds to the trajectory arrows in Fig. 3. 

Equations (31)-(34) describe a one-dimensional nonlin- 
ear eigenvalue problem for K and T with eigenvalue c. These 
equations must be solved numerically. We give the results of 
such computations in a following subsection. But first we 
prove some results concerning the direction of front propa- 
gation. 



5.1 Equations in a semi- local model 

In order to best study a nonlinear front solutions we con- 
centrate on a small finite region at a certain radius Ro in 
the disk. We assume that any kind of front structure will be 
much shorter than the local radius. We hence drop the term 
arising from the cylindrical geometry in Eqs (15)-(16). We 



5.2 Direction of front propagation 

The hypothetical front solution may not be static. But in 
which direction will it move and at what speed? If it connects 
a quiescent state to a turbulent state, will it move into the 
turbulent state (c < 0) or will it move into the quiet state 
(c > 0)? Later, it is shown numerically that either case is 
possible. In this section we try and provide some arguments 
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Figure 5. Simple cartoon of the temperature profile of two fronts. 
The red line overlaying them indicates the critical temperature 
separating the two basins of attraction of the hot/turbulent and 
cold/quiescent states. 



Figure 4. The combined heating and cooling rate as a function of 
temperature for the slaved model. Arrows indicate the direction 
in which gas will tend to move if subject to the thermal dynamics. 
We can easily see the two basins of attractions of Tq = 1 and Ta , 
with the intermediate state acting as the boundary between them. 



for why this is so and offer some analytic predictions for 
what sign c will take. 

The main ideas can be best illustrated using the simple 
'slaved' system introduced in Section 3.4. Upon imposing 
the assumptions of the previous subsection, the governing 
equation (19) becomes 



d 2 T dT 
~d!? +C ~dJ 



• r-A = o, 



Multiplying by dT/dz and integrating, one gets 
fZ (T-A)d i Td^ 



and so 



/rjaeTpde 



f T A 

sgn(c) = sgn / (r - A) dT. 



(36) 



(37) 



(38) 



This tells us that the sign of c depends on the 'area' (in T- 
space) that lies beneath the combined heating and cooling 
rates. 

The combined heating and cooling (T — A) is plotted 
in Fig. 4. Quite generally there will be two opposing contri- 
butions to the area under this curve (the integral in (38)). 
There will be a positive contribution from the portion of the 
front lying near Ta, in the active state's basin of attraction. 
And there will be a negative contribution from the portion 
of the front near the lower quiet state, in its basin of attrac- 
tion. Because the slaved model is only ID, these two basins 
are easy to see. Gas that is a little hotter than the interme- 
diate state is drawn to the upper turbulent state, while gas a 
little colder is drawn to the lower state: this is represented in 
Fig. 4 by the horizontal arrows. For the parameters chosen 
here the contribution from the active state's basin dominates 
and c > 0. 

Even if the integral in Eq. (38) were intractable, we 
can at least say that the relative sizes of the two basins of 
attraction help determine the sign of c. For example, if the 



parameter b/w is very small then Ta = (1 + w/b) 1 ^ 4 will be 
large, and consequently the basin of attraction of the upper 
active state will also be large. We would then expect the 
integral (38) to be dominated by this region and for c to be 
positive. On the other hand, if b/w is continuously decreased 
so that Ta approaches the intermediate state, then its basin 
of attraction will also decrease. At a second critical b/w 
one would expect the sign of c to change from positive to 
negative. 

These heuristic arguments can be worked through in 
some detail because (38) is integrable once we treat the heat- 
ing r as essentially a step function (cf. Section 4.1). We find 

i-T A 

(r - A)dT « i bT% - (b + w)T A +wT M m + f b, 

in which Ta ~ (1 + w/b) 1 ' 4 . This equation permits us to 
calculate the critical value of b/w at which point c = 0, as 
a function of Tmri- The equation corresponding to c = is 

(1 + w/b) 5/4 - |(w/6) Tmri -1 = 0, (39) 

which must be solved numerically. However, because b/w 
is small, a reasonable approximation can be obtained by 
expanding the first bracketed term. This gives 

(b/w) CT « (5T MR i/4)- 4 a 0.0025. 

It is easy to see that the system admits two families of MRI- 
fronts. When (b/w) < 0.0025 we have an outward moving 
front, c> 0. When 0.0042 > (b/w) > 0.0025 we have an in- 
ward moving front, c < 0. And when (b/w) > 0.0042, there 
is no front at all (because there is no active state available, 
cf. Section 4.1). Later we discuss how this parameter config- 
uration b/w depends on disk radius and, consequently, how 
an MRI-front changes as it propagates through the disk. The 
fact that there is a value for b/w for which c = means that 
there is likely a special radius at which c — 0, and hence an 
attracting radius towards which all fronts migrate. 

Finally, we may obtain a better physical understanding 
on the sign of c through the following argument. In Fig. 5 we 
plot two hypothetical MRI-front profiles for the temperature 
T. Overlying these profiles is a red line which corresponds to 
the critical temperature Ti separating the basins of attrac- 
tion of the quiet and the active states. In the first picture the 
majority of the gas lies within the active zone's basin. Qui- 
escent fluid heavily perturbed by the neighbouring front will 
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Figure 6. MRI-front profile with P r = 1, 6 = 0.001, to = 1. 
The front speed is c = 0.28Lyf2, so the front travels into the 
dead-zone, i.e. outward. 



Figure 7. MRI-front profile with P T = 1, 6 = 0.001, w = 0.457. 
The front speed is almost zero c = 9.3 X 10 — 4 Lj<£l. Note the 
extended morphology of the temperature transition in comparison 
to K. 



find itself more likely to be pushed over this critical thresh- 
old than fully turbulent fluid: the hot turbulent fluid will 
have to cool off much more than the cold fluid has to heat 
up. As this fluid heats up and becomes turbulent the front 
moves forward 'gobbling' up this gas, and hence moves fur- 
ther into the quiescent state — the large black arrow. In the 
second picture the situation is reversed and the MRI-front 
'retreats' into the turbulent gas. 



5.3 Numerical examples of MRI- fronts 

Having listed some analytic results in certain limits which 
help us understand the nature of the dead-zone/active-zone 
interface, we now compute these solutions explicitly from 
Eqs (31)-(34) using a numerical scheme. The problem is a 
4th order nonlinear eigenvalue problem in one variable £ , and 
we expect the system to be quite stiff: the profiles will not 
change much over most of the domain except near the transi- 
tion front where they will vary rapidly. As a consequence, we 
employ a relaxation method, whereby an initial guess to the 
front profile and an initial guess for c are gradually 'evolved', 
via Newton's method, to a numerical approximation of the 
correct solution (Press et al. 2007). 
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Figure 8. Front speed c as a function of the parameter w. Other 
parameters are held constant, namely 6 = 0.001. As w decreases 
so does c until at a critical point c = 0, after which c is negative 
and the front propagates inward. The critical w is approximately 
w = 0.4563. 



5.3.1 MRI-front profiles 

We show results for which b = 0.001 and w varies. MRI- 
fronts travelling either inward or outward are supported 
by the equations. Two examples are plotted in Figs 6 and 
7. In both cases the front joins the active and dead-zones 
smoothly over a characteristic lengthscale. For most out- 
ward propagating fronts, the radial scale of the transition 
region for both K and T is roughly 3Lt- However, slower 
moving fronts (in either direction) or stationary fronts ex- 
hibit a temperature transition with a long 'tail' that extends 
tens of Lt into the dead-zone. 



5.3.2 Front velocities 

Whether the front travels inward or outward its speed c is 
bounded according to \c\ < LtO. < HCl, and is thus con- 
trolled by the largest turbulent eddies; naturally, the front 
cannot travel faster than one diffusion length per orbit. The 
fastest fronts are those with larger w, i.e. originating deeper 
in the active zone at small radii. These may travel outward 
as fast as a disk scale height H per orbit: roughly 1 au in 
less than 10 years. However, c will vary as the front travels 
because of the dependence of the parameters (w especially) 
on radius, as well as the physical scales Q. and H . 



In Fig. 8, we show how the front speed c changes as we 
vary the shear parameter w, keeping b fixed. As w decreases 
so does c until at a critical point c = after which c is 
negative and the front propagates inward. At a second value 
of w the front ceases to exist: the upper turbulent state 
vanishes and the gas is no longer bistable. This behaviour 
adheres to our expectations given the simple slaved system 
of Section 5.2. There the critical value of b/w is 0.0025, at 
which point c = 0. The more general system here yields a 
value closer to 0.0022, which is in fair agreement. 

We associate smaller w values with larger radii (less 
shear) and larger w with smaller radii (greater shear). It is 
expected then that fronts originating at larger radii prop- 
agate inwards and fronts coming from smaller radii prop- 
agate outwards. In between there is a critical radius R a tt 
at which point c = and the front stays put. It follows 
that all MRI-fronts are attracted to this special radius. 
And this radius by necessity lies within the bistable region: 
R™ d < Ra.tt < i?' urb . Consider a front that begins near the 
inner edge of the bistable zone R T c ad . It will travel outward 
at some initial speed, but as it approaches R at t it will decel- 
erate; eventually it will slow to a crawl and come effectively 
to a halt near this radius. On the other hand, a front that 
begins near the outer boundary of the bistable zone _R' urb 
will propagate inwards until it too approaches i? at t , at which 
point it slows down and eventually stops. 



6 TRAVELLING FRONTS IN GLOBAL DISK 
MODELS 

In this section the physics explored in the previous sections 
are tested in simulations of the full equations (15)-(16) in a 
simple global disk model. As we have argued, an MRI-front 
usually does not linger in one place; it will by its nature 
move away to radii where the properties of the gas are going 
to be different to where it started, and consequently where 
its profile and speed will change. How this works in detail 
can only be captured by a calculation that take into account 
the radial structure of the disk. 

The geometric cylindrical terms are reinstated in (15)- 
(16) and the disk's radial structure is encoded through the 
dependence of the cooling and shear parameters b and w on 
radius R. By setting w oc Q and recognising that the cooling 
time decreases with radius, we suggest the following model 
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(40) 



where the constants wo and bo correspond to the critical 
values for which there are no active states, and at which 
point R = 7?* urb . Recall from Section 3, that this occurs 
when bo/wo = 0.0042. The power p Js we are free to set, 
but given the uncertainties in the opacity in the inner disk 
we set p — in most applications. In addition, we should 
let T eq also be a function of R because of heating from the 
central star. This aspect of the problem is neglected mainly 
for ease, and also because appreciable change in T eq should 
be limited to near the inner disk edge. 

The equations (15)-(16) are simulated with a finite- 
difference scheme and explicit first order time-stepping over 




Figure 9. Space-time diagram of the evolution of the turbulent 
amplitude it in a global model of a protoplanetary disk. The front 
propagates from near the inner edge and settles near the radius 
iiatt, which is some 65 % of ij* urb = 2 au. 



CD 



4=! 1500 



CD 
Q_ 





. 




pturb 




c 


\ Temperature 




Turbulent >v 




amplitude \. 





Radius (AU) 



Figure 10. Long time equilibrium of the temperature field in a 
global model of a protoplanetary disk. Superimposed on the figure 
is the turbulent amplitude rescaled by a factor 500 so as to be 
visible. 



a radial domain extending from 0.1 au to 5 au. Space units 
are chosen so that it' = 2 au and Lt =0.1 au. Time 
is scaled by the orbital period at _Rj; urb - Various initial con- 
ditions were trialed, including simple tanh profiles. But all 
yielded the same long term equilibrium. 

In Fig. 9 we plot a typical space-time diagram of the 
turbulent amplitude K. The front is released at a small ra- 
dius, near the inner edge, from where it propagates outward 
rapidly, covering 0.6 au in less than 5 orbits. It then slows 



down as it approaches the critical radius i? a tt and eventually 
conies effectively to a halt. In Fig. 10 we plot a snapshot of 
the long time temperature profile, obtained after 200 orbits. 
The front is completely stationary at this point. Note that, 
unlike the turbulent amplitude, the temperature T varies 
with radius even far from the location of the front. This 
is because the turbulent equilibrium state Ta is a function 
of radius through w. Thus at smaller radius the turbulent 
heating is greater, as is Ta (cf. (23)). Note that the model 
appears to overestimate T at these small radii in comparison 
with alpha models (see e.g. Terquem 2008). Superimposed 
upon the figure is the turbulent amplitude, scaled appropri- 
ately. 

In these figures, and more generally, the final location 
of the front i? att is roughly 65% of i^ urb . This could be pre- 
empted from Eq. (40), once we recognise that b/w — 0.0022 
at this radius. Then we have the estimate 



/)';. 



m 



0.0022 \ 2/(3+2p) 
0.0042 J 



(41) 



which gives ~ 0.65 when p — 0. For p > the attracting 
radius is greater. 

Finally, the front profiles we have simulated appear to 
be stable even when they slow to a static equilibrium. The 
simulations show no evidence of secondary instability and 
a subsequent dissolution into more disordered time- varying 
dynamics as witnessed in, for example, flame fronts (Zel- 
dovich and McNeill 1985). Neither is there any sign of the 
instability uncovered by Wunsch et al. (2005, 2006). If such 
instabilities exist in real protoplanetary disks, our model 
is too simple to capture them, omitting, for example, ac- 
cretion. But there is also the question of consistency. Our 
model, and the related alpha disk models, have been de- 
rived via averaging over small-scale fluctuations, and have 
thus assumed that these fluctuations are well-behaved and 
bounded. The issue of front instability and more compli- 
cated time-dependent behaviour must be settled by direct 
numerical simulations of the MRI itself. 



7 SUMMARY AND DISCUSSION 

We now summarise the contents of the paper and bring out 
its most important points. In Section 2, we made the case 
that in most protostellar disk systems, midplane gas at inner 
radii (less than about 1 au) finds itself bistable: it can be 
either quiescent or MRI-turbulent. As a consequence, the 
location of the dead/active zone boundary may be difficult 
to pin down and, moreover, the boundary profile itself not 
particularly defined. Estimates on the critical radii which 
bound this gas are not well constrained; but we expect the 
inner radius R T c lld to be (much) less than 0.1 au, and the 
outer radius i?* urb to be roughly 1 au, with its exact value 
dependent on stellar mass and accretion rate. 

In order to understand this region, a reduced model 
was introduced (Section 3) which describes the mutual cou- 
pling of the turbulent and thermal dynamics of the gas 
in the bistable region. The turbulence is modelled in a 
mean-field fashion, which assumes that its fluctuations are 
well behaved: large amplitude intermittent behaviour (which 



might be quite important) is smoothed out. The govern- 
ing equations take the form of two reaction-diffusion equa- 
tions for temperature and turbulent amplitude. After a suit- 
able rescaling, the remaining free parameters have relatively 
transparent meanings: b represents the local gas opacity, and 
w represents the local shear. 

In Section 4, we set this model to work computing ho- 
mogeneous equlibrium states, in particular showing that in 
the bistable zone of a disk there should be three equilibria: 
the dead-zone state and the active turbulent state (both sta- 
ble as expected) in addition to an intermediate state that is 
unstable and which polices the boundary between these two 
outcomes. 

In Section 5 we compute what we call 'MRI- fronts', 
which are travelling transitions connecting the active zone to 
the dead zone. We can prove a number of properties about 
these fronts. Of greatest importance is the sign of the front 
speed: we show that fronts beginning at the inner edge of the 
bistable zone propagate outward, while fronts beginning at 
the outer edge of the same zone travel inwards. There then 
exists a critical radius i? a tt at which point the MRI-front ve- 
locities are precisely zero and it is towards this radius that 
all fronts must go and finally anchor themselves. This is the 
true effective boundary between active and dead zones. Nu- 
merical simulations of the equations in a global model of a 
protoplanetary disk are given in Section 6 and these confirm 
this expectation. They also show that the front solutions are 
stable to secondary perturbations, at least within the con- 
fines of our model assumptions. 

So where does this leave us? Our simple system predicts 
that the interface between MRI-active and inactive regions 
is in fact a robust and coherent feature that relaxes into 
a steady profile. Large deviations from this steady equilib- 
rium will decay and the MRI-front will ultimately return to 
its halting radius i? a tt- This critical radius can be somewhat 
less than existing estimates, but the idea that there exists 
a well-defined boundary still applies. The picture that we 
have developed is certainly favourable to planet formation 
scenarios that rely on a stable and quasi-static dead/active 
zone interface. If the interface is as robust as suggested by 
the model, its associated pressure inhomogeneity may be 
equally robust and long-lived: precisely what is required for 
the accumulation and eventual agglomeration of solid mate- 
rial. 

However, the model has neglected (in fact, 'averaged 
away') the highly intermittent turbulent fluctuations that 
we expect such an MRI-front to suffer. The smooth profiles 
that we have generated (Figs 6 and 7) represent an aver- 
age; in reality, an MRI front will most likely be a disordered 
structure subject to potentially large amplitude perturba- 
tions. Moreover, it is not quite clear what lengthscales this 
disorder encompasses - it will probably be on a scales of 
order the transition scale itself - some fraction of H. If the 
transition region is indeed volatile and fairly broad then it 
may inhibit the development of the pressure inhomogene- 
ity required by planet formation theories. Moreover, addi- 
tional physics (accretion, for example) and realistic vertical 
structure may stimulate instability and associated secondary 
fluctuations. A sensible way to check some of these effects 
is to better describe the turbulent dynamics and move on 



from the diffusive approximation. MHD simulations of the 
MRI with a temperature dependent resistivity are needed. 
These might be undertaken in cylindrical geometry, which 
is sufficient to capture the key global effects (with shear rate 
declining with radius), while allowing acceptable resolution, 
given modern computational constraints. Such simulations 
are planned and will be presented in later work. 
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APPENDIX A: DERIVATION OF THE 
MEAN-FIELD MODEL 

Though the mean-field model used throughout the main pa- 
per can be motivated by physical heuristic arguments, its 
general lineaments can formally be derived through a sepa- 
ration of scales analysis. In this appendix we briefly sketch 
out the procedure. 



Al Governing equations 

We begin with the full equations of resistive and viscous 
MHD 



D t p = -pV ■ u, 



(Al) 



B' 



1 



Dm = -V$ - - [P + — }+ — B VB + - v n, 

p \ Stv J Airp p 

(A2) 



D t B = BVu-BVu + j?V B, 

E D t S = p(V>dis + V«t) - V • F, 



(A3) 
(A4) 



where p is volumetric mass density, u is velocity, $ is the 
gravitational potential of the star, P is pressure, B is mag- 
netic field, II is the (molecular) viscous stress tensor, r/ is the 
resistivity, E is the internal energy density, S = InPp ' is 
the entropy function and 7 the ratio of specific heats, p^dis 
comprises viscous and resistive heating, and pVcxt represents 
external heating from the star. Finally, the radiative flux is 
denoted by F. We also have denoted the total derivative by 
D t = dt + u ■ V. In writing the above we have assumed an 
ideal gas equation 



P=^pT 



(A5) 



where 1Z is the gas constant and p is mean molecular weight. 
It follows that E — P/(y — 1). We set the resistivity to be a 
function of temperature r\ = r)(T) via an appropriate form 
of Saha's equation (see Stone et al. 2000 and Balbus 2011). 
This means that we are considering thermal ionisation ex- 
plicitly. We neglect non-thermal ionisation. 



A2 Fluctuations 

We assume that equations (A1)-(A4) admit a laminar equi- 
librium state characterised by (near) Keplerian rotation: 

u = u {R) = Rn(R)e 4 , 

where R is cylindrical radius (the tiny radial drift arising 
from molecular viscous torques we omit). There will also be 
some radial and vertical structure in the variables p and P, 
which won't be listed. Finally, we assume that there is a 
weak magnetic field Bo that does not influence the equilib- 
rium balances significantly but which plays the defining role 
in the disk's stability by catalysing the MRI. 

In a given region of this disk the MRI will occur if the 
local resistivity is sufficiently low. We next assume that in 
such an unstable region turbulence will ensue, giving rise 
to both magnetic and velocity fluctuations which we denote 
by v and b. Because MRI turbulence is mainly subsonic, 
associated density fluctuations are considered to be small 
and are neglected. 

The velocity fluctuations are small compared to the Ke- 
plerian motion |v| <C \R£l\ and inhabit a well-defined range 
of lengthscales, the upper bound of which we denote by I, 
understood as either the main energy injection scale of the 
MRI, or the scale of the largest turbulent eddies. We then 
posit the following hierarchy I <C H <C R, where H is the 
disk semithickness. The first scaling is marginally true at 
best, and in fact we can only probably say I < H , but the 
stronger form is necessary for the mathematical argument 
to go through. This motivates the introduction of the small 
parameter e = l/H. 



A3 Fluctuation equations 

The magnitude of the fluctuations is defined through 

fc=I[pv 2 + b 2 /(87r)] 

which we term the 'turbulent amplitude' or 'intensity'. One 
can derive an evolution equation for k from Eqs (A1)-(A4). 



This can take multiple forms (see for example Balbus and 
Hawley 1998). We use the following: 



dtk + \7 -(vfe) 



(A6) 



where the source term q is rather ungainly. It can be written 

as 



q = v • (v • n) - v • VP 

1 

+ — 

47T 

+b • (b • Vu ) + (b + Bo) • V(b ■ v 



V • (u k) 
rjb- V 2 b + v- (v Vu ) 



v-V(Bo-b)]. 

(A7) 



This q term is obviously an exceedingly complicated function 
of k. It both controls the emergence of the MRI and its 
nonlinear saturation through turbulence. Note that it also 
depends explicitly on temperature through 77, and to a lesser 
extent P. 

In addition, we can derive an equation for fluctuations 
in internal energy or temperature. It is more convenient, 
however, to deal with the total temperature T. It is 



d t T + V -(vT) = 
(7-1) 



+ 



-R.jp 



-V-(u T) + (2-7)TV-u 



?Adis + V«xt — v 
p 



(A8) 



In both Eqs (A6) and (A8) the left hand flux term comprises 
only the advection from the turbulent velocity field. 



A4 Small scales and averages 

We wish to examine the dynamics on scales longer than the 
fluctuations, though shorter than the radial scale of the en- 
tire disk. Basically, we want to look at stuff happening in 
the bistable region which has a radial size on an interme- 
diate scale < 1 au. We hence introduce a spatial average 
over (a) azimuth, (b) the vertical extent of the disk, and 
(c) over the short radial scales of order I. This procedure 
effectively smooths out turbulent fluctuations and reduces 
them eventually to a turbulent flux, a violation of some of 
the complexities of the flow, but a procedure that furnishes 
relatively simple evolution equations for the two quantities 
of interest k and T. Note that, after averaging, these quan- 
tities depend only on t and R, where the latter coordinate 
describes radial scales longer than the turbulent eddy size I. 
In order to define proper averages over the turbulent 
fluctuation scale I at a given point in the disk, a short ra- 
dial coordinate is necessary. This we denote by x and define 
through 

R — Ro + ex 

where Ro is the radius in the disk in which we are interested. 
We may now regard any fluctuating quantity (such as v and 
b) as depending on both x and Ro (in addition to <f>, Z, 
and t). The R derivatives of these functions must then be 
replaced by Ro and x derivatives: 



_d_ 
dR 



d 

dRo 



ld_ 

t dx 



A multiple spatial average is introduced. For a quantity 



/ = f(x, i?o, 4>y Z, t) we have 

r 2n 



(f) 



1 



2 3 nHL 



dd 



H ,-L 

dZ I dxf, 

H J-L 



(A9) 



where L is an intermediate scale of order H upon which 
we wish to investigate the thermal/turbulent dynamics. It 
follows then that the averaged quantity (/) depends on only 
Ro and t. 



A5 Mean field equations 

We first attack the evolution equation for k, i.e. Eq.(A6), 
using the average defined in the previous subsection. The 
most problematic term is the flux (V • (vfc)). We first note 
the following 

r 2ir 



I d<t,{v4, k)d(f> — [vj, k]'^ — 
Jo 



d z (v z m)dz — [v z m]_ H — 0, 



where the first equality holds due to periodicity in <f>, and 
the second holds providing that velocity fluctuations go to 
zero at the disk's upper and lower surfaces. This leaves the 
R derivative term, 

d R {Rv R k)dx - — — / Rv R kdx + e~ 1 [Rv R k]^. L . 
-L 9Ro J_ L 

Recognising that I <C L, the velocity fluctuations v R arc 
uncorrelated at x — ±L. Also L < J? so that R « Ro. 
This means that the second 'surface' term is merely a source 
of 'white noise' to leading order and thus contributes no 
systematic bias over intermediate and long times. Formally 
it can be averaged out by an additional temporal average 
over short times. Here we just let it equal zero. 
In summary, 

( v -( vfe ))«-^^r(^(^ ft ))> ( A1 °) 

it oR 

and so the averaged turbulent amplitude equation becomes 
d 



d_ _L_ 

dV ' + R dR 



(Ro(v R k)) = 



(All) 



where Q = (q). This is an evolution equation for the aver- 
aged turbulent amplitude (k), which we now relabel K. 

An analogous procedure gives an equation for the aver- 
age temperature, 

d 



dV ' Ro dRo 



[Ro ({v R T) + G)} 



A, 



where the radial radiative flux O is given by 

(7-1) 



e 



( P )TZ/fi 
the 'turbulent heating' function is 

(7-1) 



(F R ), 



r = 



n/u, 



{^di s >-(2- 7 )(TV-v) 



(A12) 



(A13) 



(A14) 



and 



(7-1) 



(V»e 



+ 



2 3 HLty 



[F, 



iii 



H /(p}d<pdx , 
(A15) 



is a 'cooling' function comprising the external source and 
radiative losses from the disk surfaces. Unfortunately, in or- 
der, to get the radiative flux and cooling term in this form, it 
is necessary to treat the 1/p factor of V • F in (A8) as (1/p), 
i.e. to use a 'pre- averaged' value for the density, instead of 
properly taking account of its variation in the full integral. 
It is doubtful, however, that this introduces any apprecia- 
ble change to the cooling function A itself. Henceforth we 
suppress the angle brackets on (T). 

In order to make any progress with Equations (All) 
and (A12) we need to know how the functions Q, F, and 
A, and the fluxes O, (v R T), and (v R k) behave as K and 
T vary. Until we know this, we cannot proceed. This is the 
closure problem. We prescribe the form of these functions 
from simple but physically plausible assumptions. In Section 
3.1 some of these choices are discussed, we expand on this 
in the following subsections. 



A6 Fluxes: the diffusion approximation 

Because the disk is optically thick at the radii in which we 
are interested (~ 1 au) we are permitted to employ the dif- 
fusion approximation for the radiative flux of energy. That 
is to say that energy radiates preferentially down a temper- 
ature gradient. Consequently, 



9 



-D? d dT 



dRo' 



(A16) 



where the diffusion coefficient is a function of temperature. 
The turbulent fluxes are naturally more difficult. We 
employ the simplest turbulent closure in this instance — 
that of the eddy viscosity model — in which turbulence sim- 
ply transports gas quantities down gradients in those quan- 
tities. If both T and k were proper passive scalars then this 
assumption could be made more strongly using an additional 
multiscale method (see Frisch 1989), or a Lagrangian type 
argument (see Balbus 2011). Unfortunately they are active 
scalars and interact with and vary along the eddies that 
convect them. We believe however that this model suffices 
to get a handle on the basic qualitative behaviour in which 
we're interested. We take the following: 



(v R T)^-Dr h ^, (v R k) 
oRo 



-D f 



dK 
dRo' 



(A17) 



where the diffusion coefficents D^ Th and Dk are functions 
of the local (averaged) turbulent magnitude K. Theoretical 
issues with this approach are raised and discussed in Section 
3.5. 

Lastly we examine the relative magnitudes of turbulent 
and radiative diffusion of heat. Though MRI turbulence will 
certainly mix heat and will spread its turbulent motions, 
these transport processes have not been studied in any depth 
(see, however Hirose and Turner 2011 for an examination of 
the optically thin upper layers of a disk) . Numerical work has 
focused mainly on the transport of angular momentum and 
small particles. For instance, Carballido et al. (2005) mea- 
sure the diffusion of a passive scalar in MRI turbulence sim- 
ulations and find the diffusion coefficient to be ~ 10~ 2 H 2 Q. 
This quantity, however, depends on the magnetic field ge- 
ometry and strength and can take larger values (Johansen 



et al. 2006). Adopting a minimum mass solar nebula, we 
estimate the turbulent diffusion coefficients to be 



D f 



D\ 



ctthH Q. ~ QthlO 1 cm s , 



at 1 au, where ath is the 'turbulent alpha' for thermal trans- 
port < 10~ 2 . These coefficients drop to zero, of course, in 
the quiescent state. 

Radiative diffusion, on the other hand, is easier to con- 
strain. We have 



D 



■ad 16(7 - 1)o- S bT 4 



3K,pP 



2.39 x 10 1 



/ T \ 3 /cmV 



ViooKy I K 



in which we have set p ~ 10 g cm" . The Stefan- 
Boltzmann constant is ctsb- It follows that in a hot turbu- 
lent state (~ 1000 K), radiative diffusion is comparable to 
turbulent diffusion. In the quiescent state, however, radia- 
tive diffusion dominates, as then D^ T goes to zero. Note 
that beyond the dust sublimation temperature > 1600 K 
the opacity will vary steeply and, as a consequence, radia- 
tive diffusion will certainly dominate. However, it is expected 
that turbulent heating sustains such temperatures at very 
small radii (< 0.1 au; see Terquem 2008, for example). In the 
present work we do not attempt to model the complexities 
of the transition from one regime to the other and instead 
treat the diffusion coefficients as constants. 



A7 Reaction terms 

We address A first, which determines radiative equilibrium: 
the temperature balance in the absence of turbulent dissi- 
pation. Radiative losses on the disk surfaces can be approx- 
imated by 



2 3 HLn 



PZTZ P Li 

I / [F Z ] H H d<j>dx : 

J -L 



2ctsbT but 



where T sur f is a radially and azimuthally averaged effective 
surface temperature. We relate T su rf to the midplane tem- 
perature T c through T BUl{ = 4T 4 /(3t c ), as is customary, 
where r c is half the total optical thickness. The central tem- 
perature T c is then assumed to be proportional to the verti- 
cally averaged temperature T. Putting this into (A15) and 
rearranging gives 



A = 6(T 4 -T 4 q 



(A18) 



where T cq is the temperature at radiative equilibrium and 
6 is some free parameter associated with the opacity of the 
gas. It controls, in basic terms, the speed at which radiative 
equilibrium is enforced. 

The MRI terms are represented by the Landau operator 
Q = sK — aK 2 , where s is the linear growth rate of the MRI. 
It is assumed that temperature variations only intrude on 
the model via this linear forcing term. Smaller temperatures 
mean greater resistivities, and hence lower (or zero) growth 
rates. From consideration of the linear dispersion relation of 
the resistive MRI, Lesaffre et al. (2009) accounted for this 
via s = sq[1 — v{T)] where so is the ideal MHD growth 



rate, r\ is the non-dimensionalised resistivity. From the Saha 
equation a plausible form for r\ might be 

770c T~ 1/4 cxp(T7T), (A19) 

for T* a constant. However, the extreme gradients, and mas- 
sive decay rates that can ensue, mean that the prescription 
(A19) has terrible numerical properties. A more convenient 
model for r\ is 



77 oc 1 - tanh[6(T - T M Ri)/T e( 



(A20) 



at least in the context of our reduced system. The tanh func- 
tion captures the 'switch'-like behaviour of the Saha equa- 
tion at Tmri but does not exhibit a negative and divergent s 
when T is low. The choice of coefficient 6 gives the best com- 
parison with (A19) in the regimes of interest. For instance, 
the maximum relative error that emerges in s is some 6%. 



APPENDIX B: LINEAR STABILITY OF 
HOMOGENEOUS STATES 

We investigate purely homogeneous perturbations so that 

K = K + K 1 (t), T = T Q + T 1 (t) 

where Ko and To is the equilibrium in question and K\ and 
T\ are small perturbations. The linearised equations of (15) 
and (16) are 

<Ttfi=«o#i-*o#oTi-2tfotfi, (Bl) 

0-T1 = wKi -46T 3 Ti, (B2) 

in which we have decomposed T\ and K\ into Fourier modes 
oc e at , where a is a linear growth rate (not to be confused 
with the MRI growth rate s). Also so and s correspond to s 
and ds/dT evaluated at T — Tq. The following second order 
dispersion relation ensues: 



a 2 + {Kl + 46T 3 )a + (46 T 3 - ws' )K = 0. 



(B3) 



Linear stability of the state (Ko, To) is assured if both the 
coefficient of a and the last term are positive. The first con- 
dition is always satisfied, which leaves us with the main sta- 
bility criterion: 

4foT 3 -ws'o > 0. (B4) 

This can be rewritten into a more illuminating form by using 
the cooling function A and the heating function at equilib- 
rium 

T(T) = wK = (lis. 

Then the stability criterion becomes 



dA dT 
dT > dT' 



(B5) 



which says that the cooling rate must outstrip the heating 
rate if the state is to be stable, in accord with the intuitive 
argument given in Section 4.2. 



